%% US Public Debt and Safe Asset Market Power
%% Jason Choi, Rishabh Kirpalani, and Diego Perez
%% Nov 24, 2024

%% Solve Competitive Equilibrium

%----------------------------------------------------------------
% 0. Housekeeping
%----------------------------------------------------------------

close all

%----------------------------------------------------------------
% 1. Defining variables
%----------------------------------------------------------------

// Endogenous Variables
var b rb rkstar rk krw_star kus_star krw kus kstar k wstar w c_rw c_us drdb spread vrw vus y dMrwdb;

// Exogenous Variables
var nnu oomega A Astar;

// Shocks
varexo eps_nnu eps_oomega eps_A;

// Parameters
parameters ggamma bbeta eeta llambda aalpha Astarbar Abar iiota iiota_star ddelta_rw ddelta_us
  nnu_bar oomega_bar rrho_nnu rrho_oomega ssigma_nnu ssigma_oomega rrho_A ssigma_A rrho_Astar ssigma_Astar
  b_ce rb_ce rkstar_ce rk_ce krw_star_ce kus_star_ce krw_ce kus_ce kstar_ce k_ce
  capKstar_ce capK_ce wstar_ce w_ce crw_ce cus_ce spread_ce vrw_ce vus_ce y_ce;

%----------------------------------------------------------------
% 2. Calibration
%----------------------------------------------------------------

% // Parameters
disp('% Epsilon = 2.2, Lambda = 2') 
ggamma = 2;
bbeta = 0.9886;
eeta = 0.545;
llambda = 2;
aalpha = 0.3;
Astarbar = 0.9254;
Abar = 0.8154;
iiota = 0.9070;
iiota_star = 0.7939;
ddelta_rw = 0.1;
ddelta_us = 0.1;
nnu_bar = 0.0041325;
oomega_bar = 0.0165;
rrho_nnu = 0.99;
ssigma_nnu = 0.01;
rrho_oomega = 0.95;
ssigma_oomega = 0.3;
rrho_A = 0.95;
ssigma_A = 0.01;
rrho_Astar = rrho_A;
ssigma_Astar = ssigma_A;

% Analytic Steady State (Competitive Equilibrium)
b_ce = (oomega_bar/nnu_bar)^(1/(eeta-1-llambda));
rb_ce = 1/bbeta - nnu_bar*(b_ce)^(eeta-1) - 1;
rkstar_ce = 1/bbeta	+ ddelta_rw - 1;
rk_ce = 1/bbeta	+ ddelta_us - 1;
krw_star_ce = ((aalpha*(1-iiota_star)*Astarbar*((1/bbeta+ddelta_rw-1)/(aalpha*iiota_star*Astarbar))^((aalpha*(1-iiota_star)-1)/(aalpha*(1-iiota_star))))/(1/bbeta+ddelta_us-1))^((aalpha*(1-iiota_star))/(1-aalpha));
kus_star_ce = ((1/bbeta+ddelta_rw-1)/(aalpha*iiota_star*Astarbar))^(1/(aalpha*(1-iiota_star)))*krw_star_ce^((1-iiota_star*aalpha)/(aalpha*(1-iiota_star)));
krw_ce = ((aalpha*iiota*Abar*((1/bbeta+ddelta_us-1)/(aalpha*(1-iiota)*Abar))^((aalpha*iiota-1)/(aalpha*iiota)))/(1/bbeta+ddelta_rw-1))^((aalpha*iiota)/(1-aalpha));
kus_ce = ((1/bbeta+ddelta_us-1)/(aalpha*(1-iiota)*Abar))^(1/(aalpha*iiota))*krw_ce^((1-(1-iiota)*aalpha)/(aalpha*iiota));
kstar_ce = krw_star_ce + krw_ce;
k_ce = kus_star_ce + kus_ce;
capKstar_ce = krw_star_ce^iiota_star*kus_star_ce^(1-iiota_star);
capK_ce = krw_ce^(1-iiota)*kus_ce^iiota;
wstar_ce = Astarbar*(1-aalpha)*(capKstar_ce)^aalpha;
w_ce = Abar*(1-aalpha)*(capK_ce)^aalpha;
crw_ce = wstar_ce + (rkstar_ce-ddelta_rw)*kstar_ce + nnu_bar/eeta*(b_ce)^eeta + rb_ce*b_ce;
cus_ce = w_ce + (rk_ce-ddelta_us)*k_ce - oomega_bar/(1+llambda)*(b_ce)^(1+llambda) - rb_ce*b_ce;
drdb_ce = 0;
dMrwdb_ce = 0;
nnu_ce = nnu_bar;
oomega_ce = oomega_bar;
A_ce = Abar;
Astar_ce = Astarbar;
spread_ce = (rk_ce-ddelta_us-rb_ce);
vrw_ce = crw_ce^(1-ggamma)/(1-ggamma)/(1-bbeta);
vus_ce = cus_ce^(1-ggamma)/(1-ggamma)/(1-bbeta);
y_ce = A_ce*(kus_ce^iiota*krw_ce^(1-iiota))^aalpha;

%----------------------------------------------------------------
% 3. Model
%----------------------------------------------------------------

model;

c_rw^(-ggamma) = bbeta*c_rw(+1)^(-ggamma)*(nnu*b^(eeta-1)+1+rb);
c_rw^(-ggamma) = bbeta*c_rw(+1)^(-ggamma)*(1-ddelta_rw+rkstar);
c_rw + kstar + b = wstar + (1-ddelta_rw+rkstar(-1))*kstar(-1) + nnu(-1)/eeta*(b(-1))^eeta + (1+rb(-1))*b(-1);

c_us^(-ggamma) = bbeta*c_us(+1)^(-ggamma)*(oomega*b^llambda+1+rb+drdb*b);
c_us^(-ggamma) = bbeta*c_us(+1)^(-ggamma)*(1-ddelta_us+rk);
c_us + k - b = w + (1-ddelta_us+rk(-1))*k(-1) - oomega(-1)/(1+llambda)*(b(-1))^(1+llambda) - (1+rb(-1))*b(-1);

rk = Astar*aalpha*(1-iiota_star)*krw_star^(aalpha*iiota_star)*kus_star^(aalpha*(1-iiota_star)-1);
rkstar = Astar*aalpha*iiota_star*krw_star^(aalpha*iiota_star-1)*kus_star^(aalpha*(1-iiota_star));
rk = A*aalpha*iiota*krw^(aalpha*(1-iiota))*kus^(aalpha*iiota-1);
rkstar = A*aalpha*(1-iiota)*krw^(aalpha*(1-iiota)-1)*kus^(aalpha*iiota);
wstar = Astar*(1-aalpha)*krw_star^(aalpha*iiota_star)*kus_star^(aalpha*(1-iiota_star));
w = A*(1-aalpha)*krw^(aalpha*(1-iiota))*kus^(aalpha*iiota);

drdb = 0;
dMrwdb = 0;

k = kus + kus_star;
kstar = krw + krw_star;

log(nnu) = (1-rrho_nnu)*log(nnu_bar) + rrho_nnu*log(nnu(-1)) + ssigma_nnu*eps_nnu;
log(oomega) = (1-rrho_oomega)*log(oomega_bar) + rrho_oomega*log(oomega(-1)) + ssigma_oomega*eps_oomega;
log(A) = (1-rrho_A)*log(Abar) + rrho_A*log(A(-1)) + ssigma_A*eps_A;
log(Astar) = (1-rrho_Astar)*log(Astarbar) + rrho_Astar*log(Astar(-1)) + ssigma_Astar*eps_A;

spread = (rk-ddelta_us-rb);

vrw = c_rw^(1-ggamma)/(1-ggamma) + bbeta*vrw(+1);
vus = c_us^(1-ggamma)/(1-ggamma) + bbeta*vus(+1);

y = A*(kus^iiota*krw^(1-iiota))^aalpha;

end;

%----------------------------------------------------------------
% 4. Computation
%----------------------------------------------------------------

initval;
  b = b_ce;
  rb = rb_ce;
  rkstar = rkstar_ce;
  rk = rk_ce;
  krw_star = krw_star_ce;
  kus_star = kus_star_ce;
  krw = krw_ce;
  kus = kus_ce;
  kstar = kstar_ce;
  k = k_ce;
  wstar = wstar_ce;
  w = w_ce;
  c_rw = crw_ce;
  c_us = cus_ce;
  drdb = drdb_ce;
  dMrwdb = dMrwdb_ce;
  nnu = nnu_ce;
  oomega = oomega_ce;
  spread = spread_ce;
  A = A_ce;
  Astar = Astar_ce;
  vrw = vrw_ce;
  vus = vus_ce;
  y = y_ce;
end;

resid;
check;

shocks;
  var eps_nnu = 1;
  var eps_oomega = 1;
  var eps_A = 1;
end;

set_dynare_seed('default');
stoch_simul(order=2,nograph,noprint,periods=100000,drop=1000,pruning);

%----------------------------------------------------------------
% 5. Generate moments
%----------------------------------------------------------------

spread_path = (rk-ddelta_us-rb)*100;
var_sp = var(spread_path);
auto_sp = autocorr(spread_path);
cost = oomega./(1+llambda).*(b).^(1+llambda);
benefit = (eeta-1).*b.*(log(b)-nnu);
var_by = var(b./y);
auto_by = autocorr(b./y);
corr_pq_by = corr(spread_path,b./y);

moments = [mean(b./y) mean(spread_path) var_by var_sp corr_pq_by auto_by(2) auto_sp(2)]';
data_mom = [0.41 0.62 0.03 0.086 -0.56 0.96 0.70]';
rowNames = {'Mean b/y','Mean sp','Var b/y','Var sp','Corr (b/y,sp)','Autocorr b/y','Autocorr sp'};
colNames = {'Model Moments','Data Moments'};
TableA0 = array2table([moments data_mom],'RowNames',rowNames,'VariableNames',colNames)

%----------------------------------------------------------------

b_ce_lh = mean(b);
rb_ce_lh = mean(rb);
spread_ce_lh = mean(spread);

oo_ce_lh = oo_;
M_ce_lh = M_;
options_ce_lh = options_;

save ce_lh_save oo_ce_lh M_ce_lh options_ce_lh b_ce_lh spread_ce_lh rb_ce_lh;
